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Abstract 

The compact form of the discontinuous Galerkin 
method allows for a detailed local analysis of the 
method in the neighborhood of the shock for a non- 
linear model problem. Insight gained from the anal- 
ysis leads to new flux formulas that are stable and 
that preserve the compactness of the method. Al- 
though developed for a model equation, the flux 
formulas are applicable to systems such as the Eu- 
ler equations. This article presents the analysis for 
methods with a degree up to 5. The analysis is ac- 
companied by supporting numerical experiments us- 
ing Burgers’ equation and the Euler equations. 

Introduction 

The discontinuous Galerkin method is being de- 
veloped as a means for obtaining a high-order shock 
capturing-capability on unstructured meshes. This 
capability is an important step toward achieving an 
efficient and robust method for aeroacoustic appli- 
cations. The objective of this work is to determine 
whether a local eigenvalue analysis can predict the 
instabilities that have been previously observed, and 
to use insight gained from the analysis to suggests 
ways to eliminate the instability. 

In reference 1, the discontinuous Galerkin method 
was formulated in a quadrature-free form that re- 
duced both the computational effort and storage re- 
quirements. Reference 2 described the implemen- 
tation of boundary conditions, including the treat- 
ment of curved walls and nonreflecting boundary 
conditions. In these works, the method was de- 
scribed in detail, and numerical results for scalar 
advection and for the Euler equations were shown 
to demonstrate the accuracy and robustness of the 
method. These studies showed numerically^ that 
piecewise linear and piecewise quadratic discontinu- 
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ous Galerkin methods were stable without the use of 
either limiters or added dissipation when applied to 
the nonlinear Burgers’ equation for a shocked case. 
However, higher order methods diverged immedi- 
ately after shock formation. Also in reference 1, it 
was observed in numerical test that the quadratic 
form was stable only if the flux integral term was 
evaluated exactly. 

The discontinuous Galerkin has a number of fun- 
damental properties essential to any robust shock 
capturing method. In a series of papers, Cock- 
burn and ShT and Cockburn et alf’ ® discussed 
the discontinuous Galerkin method with the use of 
approximate Riemann solvers, limiters, and total- 
variation-diminishing (TVD) Runge-Kutta time dis- 
cretizations for nonlinear hyperbolic problems. Ref- 
erence 5 presents the design of a limiter that applies 
to general triangulations, maintains a high order of 
accuracy in smooth regions, and guarantees maxi- 
mum norm stability. Jiang and Shu® also proved 
that the discontinuous Galerkin method satisfies a 
local cell entropy inequality for the square entropy 
r/(U) = IJ 2 , for arbitrary triangulations in any space 
dimension, and for any order of accuracy. This proof 
trivially implies the L 2 stability of the method for 
nonlinear shocked problems in the scalar case. 

Flux limiting has been demonstrated as a means 
for stabilizing the shocked case;®’ ^ however, this ap- 
proach tends to reduce the formal accuracy and to 
diminish the compactness of the method. This work 
presents a local eigenvalue analysis that predicts the 
stability, or instability, of the method for the Burg- 
ers’ equation when a shock is present. The analysis is 
possible only because of the compact nature of the 
discontinuous Galerkin method. The analysis also 
provides physical insight into possible causes of the 
instability and leads to flux formulas that eliminate 
the instability. 

The first section describes the discontinuous 
Galerkin method applied to the nonlinear Burgers’ 
equation and provides the rationale for the analysis. 
The second section describes in detail the analysis 
and the conclusions drawn from it. The third sec- 
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tion presents numerical experiments with the new 
flux function using the scalar Burgers’ equation. The 
last section presents numerical experiments with the 
Euler equations. 


Discontinuous Galerkin Method 

Application of the semidiscrete form of the dis- 
continuous Galerkin method to the scalar hyperbolic 
equation of the form 


dU 


9F(U) _ 

dx 


(0 < x < 1, t > 0) (1) 


begins by partitioning the domain into non- 
overlapping elements that cover the domain: 0 = 
x 0 < ■ ■ ■ < Xi < Xi - |_i < ... < 1. Within each el- 
ement, the solution is approximated by a subspace 
that is defined local to the element. In the present 
work, the subspace will be the set of polynomials 
of degree < n: S = {l,^,^ 2 ,...,^ n }. Within each 
element, x = aq + A Xi A = (*i-i + aq)/ 2, 

and Aaq = aq — aq_i. The solution in element i is 
approximated by 


u (x, t) & Vi(t,t) = 

j = o 


for Xi- 1 < X < Xi. 

The evolution of the dependent variables Vij(t) 
is governed by the projection of equation (1) onto 
the same subspace that is used to approximate the 
solution: 



dF(Vj) 

dx 


d{ = 0 


(0 < j < n) 


Integration by parts yields the weak conservation 
form of the equation, which is used in the numer- 
ical solution: 


' dt 4 


• Aa 


2 ii j ~ 1 F{v i ) di 

i 


+ (hYF(xi)-(-^F(x l - 1 ) =0 ( 2 ) 


and the second integral is evaluated exactly. Formal 
error convergence properties are obtained in smooth 
regions with the flux truncated to a polynomial of 
degree n + 1. For nonlinear problems with shocks, 
stable solutions of the quadratic case ( n = 2) were 
obtained only when the flux was evaluated exactly. 

At the element boundaries, where the solution is 
not unique, the F terms are evaluated based on the 
solution that is “upwind” (in the sense of charac- 
teristic directions) of the element boundary. F is 
generally written in the form of a Riemann flux as 
F(xi ) = F r (Ul,Ur), where Ul and Ur denote so- 
lutions from the left and right sides of the element 
boundary, respectively. Hence, adjacent elements 
communicate with one another only through F. As 
a consequence, each element can be thought of as 
a somewhat autonomous entity that responds to in- 
flow boundary conditions provided by the element (s) 
“upwind” of itself. 

Elements Near a Shock 

Under somewhat idealized conditions, the stabil- 
ity of the overall solution depends heavily, if not en- 
tirely, on the behavior of the method in the element 
that contains the shock. To show this, consider the 
model problem in which F(U) = ^U 2 and equation 
(1) is supplemented with the following initial and 
boundary conditions: 

17(0, t) = -U(l,t) = 1 


U{x, 0) 


1 x < X s 
-1 x > X s 


The approximate Riemann flux commonly used with 
this model problem is of the form 

F r {Ul, Ur) = [F(Ur) + F(U l ) - A (U R - U L )] /2 

where A > \Ur + Ul|/ 2. If A = \Ur + Ul|/ 2, then 
this flux is an exact splitting 


F r (Ul,Ur) 


F(U l ) if A > 0 
F{U r ) if A < 0 


or 

Fj-vm 

+ Ax- 1 [(i) 2 F(*0 - (-i)-^-!)] = 0 (3) 
for j = 0, 1, 2, . . . , n. 

The Cj ■ U'(0 term, which contains the two in- 
tegral terms of equation (2), depends only on the 
solution within the element. In the quadrature- free 
forint the flux is approximated as a polynomial in 


and also has the following property commonly asso- 
ciated with Roe’s flux : 4 

A (U R -U l ) = ± [F(U r ) - F(U l )\ . 

Next, assume that an exact discrete solution exists 
and that the solution is sufficiently accurate such 
that the sign of A for the discrete solution is the 
same as the sign of the exact solution (i.e. , sign(U) = 
sign(A)). Also let I denote the index of the element 
within which the shock lies (i.e., */_i < X s < */). 
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Then for all element boundaries to the left of element 
/, A > 0 and F(x{) = -F(Vj(±)). Similarly, for all 
element boundaries to the right of element I , A < 
0 and F(xi) = F(Vi+i(— f)). Hence, equation (3) 
becomes 

C r v t + Ax~ 1 (±) j F(Vi(±)) 

= Ax^l-^y F(Vi-i(±)) (*•</)( 4) 

c r v t = -Ax - 1 [(i)^(v;- + 1 (-i)) 

-{-h) j F W-i(*))] (i = /) (5) 

£j ■ Vi - A x-y^Fivy-i)) 

= -Ax-y-L y F ( V l+1 (-^) (*•>/) 

( 6 ) 

When considered separately, the subdomains to 
the left and to the right of element /, which are 
governed by equations (4) and (6), respectively, are 
stable and well behaved (and trivially satisfied by 
the initial conditions). Hence, all elements to the 
left and to the right of element I can be thought of 
as simply supplying boundary conditions to element 
F Thus, the stability of the overall problem is con- 
tingent on the existence and stability of the solution 
within the element that contains the shock. 


Analysis of Shock Containing Element 

Motivated by the previous discussion, we consider 
the behavior of equation (5), in which the approx- 
imate Riemann flux terms F are evaluated at the 
state specified by the boundary conditions (the ele- 
ment subscript has been dropped for clarity): 

fi .dV i 

I J^dt-Ax - 1 I 

~ 2 ~2 

= -Ax- 1 [( i yF(-i)-(- i yF(i)] 


Linearization of the integral flux term about a spec- 
ified test solution Vo gives 

<9V 

M— - Aa; _1 AV = -AaT^S + Fo - AV 0 ) (7) 

where* M = [m jtk ], V = [vj], A = f£, F = [fj], 
S = [ 8j ], 
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0 for j even 
A for j odd 


Matrix and vector indices range from 0 through n. 


and 

1 _ 

fj = J\jZ j - 1 F(v(Z))dt 

Equation (7) is stable if the eigenvalues of M _1 A 
lie in the left-hand plane, with one exception. Note 
that the first row of A is zero, which produces at 
least one zero eigenvalue. However, the equation for 
j = 0 is simply the constraint that the time rate of 
change of the average value of the solution depends 
only on the flux imbalance, which in this analysis is 
specified. Thus, this particular eigenvalue does not 
affect the stability of the system. To avoid confusion, 
we examine the eigenvalues of a modified matrix A 
in which this zero eigenvalue has been eliminated. 
The modified equation is obtained by discarding the 
j = 0 equation and eliminating the vo component 
of the solution from the remaining equations by re- 
quiring that the average value of V be a specified 
constant. 

Because equation (7) is nonlinear, the eigenval- 
ues cannot be evaluated until the solution is known. 
However, the exact solution is not known and may 
not even exist. For the purposes of this analysis, the 
eigenvalues are evaluated by using a test solution 
Vo that is obtained by projecting the initial solution 
onto the solution subspace S. This test solution is 
given by 



* ( j d(- 2 ( j d( (8) 

I J V 


where = (X s — x) /Ax is the location of the shock 
within the element in terms of element coordinates. 

The A and its eigenvalues cr are readily computed 
for any shock position within the element. Fig- 
ures 1(a)- 1(e) show the real part of the eigenval- 
ues Re(cr) as a function of shock position for n = 
1 through 5, respectively. In all cases, all eigenval- 
ues are zero when the shock is on the boundary. For 
n < 2, all eigenvalues have a negative real part for 
— ■§■ < <i. For n > 2, the real part of the eigen- 

values becomes positive when the shock is close to 
the boundary of the element (£ s Pd ±f); and hence, 
the method is unstable for these cases. This result 
is consistent with the numerical results of reference 
1. 

When V is a polynomial of degree n, then F(V) is 
a polynomial of degree 2 n; however, the flux polyno- 
mial can be truncated to degree n + 1 with no formal 
loss of accuracy. In reference 1 , it was observed that 
the solution obtained with n = 2 was stable when the 
flux polynomial was evaluated exactly, and unstable 
when the flux polynomial was truncated. (Note that 
truncation of the flux is not an issue in the n = 1 
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case because n + 1 = 2 n.) The eigenvalues shown 
in figures 1(a)- 1(e) were computed using the exact 
flux (no truncation). Figure 2 shows the eigenvalues 
for n = 2 in which the degree of the flux polyno- 
mial is truncated to n + 1 = 3. In this case, the 
eigenvalues are positive for shock positions near the 
boundary. Again, the analysis agrees with the nu- 
merical results. 

A physical understanding of why the element may 
become unstable can be obtained by looking at the 
test solution as the shock approaches the element 
boundary. Figures 3(a) and 3(b) illustrate typical 
solutions for n = 2 and n = 3, respectively. As the 
shock approaches the right boundary, the solution 
on that boundary becomes positive. As this occurs 
the physics within the element dictate the movement 
of “mass” into the neighboring element. However, 
the approximate Riemann fluxes that are commonly 
used do not allow the flux to switch until the average 
velocity on the element boundary becomes positive, 
and this change cannot take place until an overshoot 
has occurred (i.e. , V(±) > 1). 

Alternate Flux Formulas 

If the conjecture described above is accurate, then 
the instability can perhaps be eliminated by chang- 
ing the definition of the approximate Riemann flux 
such that it becomes dependent on the solution 
within the shocked element as the shock approaches 
the edge. For the purposes of this analysis, let the 
approximate Riemann flux at an element boundary 
be defined as the sum of all contributions that leave 
the elements on either side of the element boundary: 

T r {U l ,U r ) = max(0, U L )U L /2 

+ min(0, Ur)U r /2 (9) 

The equation for the element that contains a shock 
becomes 

s .dV fi . , 

Ji—dt-Ax - 1 / 

2 J ~ 2 

+ AX- 1 [s+(iyF(V(l/2)) 

- 8 -(-iy F (v(- 1 / 2 ))] 

= -A*- 1 [(i)^(-l)-(-i)^(l)] (10) 

where = 1 if ±R(±l/2) > 0 and = 0 oth- 
erwise. When y 0, the j = 0 equation is not 
degenerate and the eigenvalues of the full matrix A 
must be considered. The eigenvalues of the A or A 
matrices produced by the above set of equations are 
shown in figures 4(a)-4(e) . Because the eigenvalues 
are unchanged when the shock position is not near 
the element boundary, these figures focus on a small 


region near £ s = 0.5. The real part of the eigenval- 
ues is positive for all shock positions for n < 4. For 
n = 5, one eigenvalue becomes positive in a small 
region near £ s = 0.487 suggesting a possible insta- 
bility. 

The flux given in equation (9) is stable for n < 5 
and has most of the properties desired of a con- 
servative flux: it is continuous, monotone, and 

T r (U,U) = F(U). However, some concern exists 
regarding the property that F r (Ul,Ur) can ex- 
ceed both F(Ul) and F(Ur) in the region Ur > 0, 
Ur < 0. For example, when Ur = —Ur = u, then 
F r (u, —u) = 2 F(u). This can be remedied either by 
limiting the flux or by taking an average of the con- 
tributions from the left and right elements (but only 
in regions were Ur ■ Ur < 0). Note that a limiter 
of the form T ^ = mm(F R (U L ,U R ), F(U l ), F(Ur)) 
does not affect the formal accuracy of the discontin- 
uous Galerkin method. The average can be either a 
simple algebraic average or a weighted average that 
depends on the relative magnitudes of the character- 
istic speeds; however, the algebraic averaging results 
in a flux that is a discontinuous function of Ur and 
Ur. 

Both limiting and weighted averaging causes the 
contribution of the flux to the shocked element to de- 
pend on the solutions in both the shocked element 
and the neighboring element. Consequently, an anal- 
ysis of just the shocked element would be of ques- 
tionable value. Such a case would require at least an 
analysis of two coupled elements; however, this anal- 
ysis is left for future work. In the case with algebraic 
averaging, the degree to which the approximate Rie- 
mann flux contributes to the solution in the shocked 
element depends only on the solution within that 
element, and the current single-element analysis is 
still possible. The equations in the shocked element 
are similar to those given in equation (10) but with 
= 1/2 if ±R(±l/2) > 0 and = 0 otherwise. 
(Note that the right-hand side also is different but 
does not enter into the analysis.) The eigenvalues 
for n = 2 and 4 are shown in figures 5(a) and 5(b), 
respectively. Both cases are stable have eigenvalues 
and are similar to the non-averaged case. Similar 
results are found for values of n < 5 while the n = 
5 case is still unstable. 

Numerical Test 

Numerical tests of the discontinuous Galerkin 
method with the modified approximate Riemann 
flux are presented in figures 6 and 7. The initial con- 
dition used for these tests is U(x, 0) = U + cos (27T*) , 
where U is either 0 or 1/2. With this initial solution, 
a shock begins to form shortly after t = 0.15. The 
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solution is evolved in time by using a third-order 
Runge-Kutta method, as described in reference 1. 

Figure 6 (a) shows the evolution of the solution up 
to t = 0.4 for n = 3 and with the domain parti- 
tioned into 20 uniformly spaced elements. Figure 
6 (b) shows solutions at t = 0.4 for n = 0 through 5. 
This grid resolution places the exact shock position 
on an element boundary. The solutions are smooth 
and stable for all values of n tested. Figure 6 (c) 
shows a similar case with 22 elements which places 
the shock at the center of an element. For n > 1, 
the solutions have considerable oscillations that are 
not aesthetically pleasing, but the oscillations are 
bounded in time and are confined to the element that 
contains the shock and the elements on either side of 
that element. Furthermore, if the element average is 
plotted, as in figure 6 (d), then the element-averaged 
solution is monotone in every case. Figures 7(a) and 
7(b) show the moving shock case that results when 
U = 1/2. Figure 7(a) shows the local solution for t 
= 0.4, n = 0 through 5, and figure 7(b) shows the lo- 
cal and element-averaged solution for n = 3. Again, 
the oscillations do not increase with n, the oscilla- 
tions are confined to the neighborhood of the shock, 
and the element-averaged solution is quite well be- 
haved. Tests that employed some form of averaging 
when Ul • Ur < 0 showed no significant change in 
the solutions. 

Euler Equations 

In this section, the quadrature-free from of discon- 
tinuous Galerkin is applied to the Euler equations 
given by 



p 


pu 

u = 

pE 

, F = 

u(pE + P) 


pu 


pu 2 + P 


where P = (7 — 1 )(pE — ±. pu 2 ). 

The first requirement for the implementation of 
the quadrature-free form of discontinuous Galerkin 
is that the flux must be written as a polynomial. But 
unlike Burgers’ equation, an exact polynomial form 
for the Euler flux does not exist. When rewritten in 
terms of the dependent variables p, pE , and pu, the 
flux involves the terms 

_ (pu)(pE) _ (pu) 2 _ (pu) 3 

J 1 i J 2 i J 3 o 

P p p 2 

Polynomial approximations of the form 

n + 1 


an obtained by projection by solving 

f l/2 z-1/2 

/ epfidz=l eiudf, 

j - 1/2 j - 1/2 

for j = 0, 1, ..., n + 1, where Ri = (pu)(pE), 
(pu ) 2 , and (pu)f 2 for i = 1, 2, and 3, respectively. 
The equations are linear in the unknowns //*,, and 
each equation set has the same matrix on the left- 
hand side of the equation. 

For initial tests, the approximate Riemann flux is 
obtained by splitting F according to the sign of the 
eigenvalues of the Jacobian of F . Thus, if 

F = AU = T[A]T _1 G 

then 

t(u l ,u r ) = (l[\+]l- 1 ) Ul u l 

+ (T [a - ] L~ 1 ) Ur IJr 
= F+(U l ) + F-(Ur) 

The test problem is periodic on the domain 0 < 
x < 1 with initial conditions 

p = 1 + 0.1 sin(27r*) 

P = P 1 , a = Vl P /P 
u = Moo ± ^rr(a - ^fl) 

pE = P/( 1 - 1) + pu 2 / 2 

that produce a single isentropic acoustic wave that 
propagates at an average speed of u Q 0 ± The 

initial discrete solution is obtained by fitting a poly- 
nomial through the exact initial solution at the 
Chebyshev-Gauss-Lobatto points within each ele- 
ment. 

The acoustic wave propagates in the indirection 
(relative to «oo), steepens to form a shock, and then 
gradually decays in time. Figures 8 (a) and 8 (b) 
show solutions at t = 3 for right-traveling waves with 
«oo = 2.0 an d 0.5^^; respectively. At these val- 
ues for Mqo, each eigenvalue of the Jacobian A is of 
a constant sign throughout the domain. Thus, the 
functional form of the flux is the same everywhere, 
and the situation in which A u L > 0 and A jj r < 0 
does not occur for any eigenvalue. As in the results 
with Burgers’ equation, the solution is stable, the 
oscillations are confined to the neighborhood of the 
shock, and they remain roughly the same magnitude 
as the formal order of the method n is increased. 

With Mqo = the shock is stationary, and 

the sign of the eigenvalue associated with the right- 
traveling wave changes across the shock. As shown 
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in figure 9, the solutions are considerably more oscil- 
latory in this case, and a distinct difference is evident 
between solutions obtained with odd and even val- 
ues of n. In the region on the left side of the shock, 
the oscillations do not grow as n is increased, but 
they are not confined to the shock region as in the 
moving-shock case. In the region on the right side of 
the shock, results with odd values of n are smooth, 
but results with even values of n are highly oscilla- 
tory. The reasons for this unusual dependence on 
n is not known at this time. Results obtained with 
this simple flux are stable but lack sufficient accu- 
racy, and further work is needed. 

The appearance of the solution can be improved 
either by adding dissipation or by applying a filter 
to the solution at each time step. Efforts to develop 
high-order dissipation terms have not been success- 
ful thus far; however, filters that are both compact 
and robust have been developed. The filter has the 
form V = T - 1 dTV, where T is the constant matrix 
that relates the coefficients of a general polynomial 
to the coefficients of an equivalent Legendre expan- 
sion. The diagonal matrix d = (eo, e\, . . . , e m , . . .) 
contains the damping coefficients. Figure 10 shows 
results with e 8 - = e l for several constant values of 
e. This filter is capable of removing oscillations; 
however, the accuracy is reduced to first order for 
smooth flows. Development of high-order Liters and 
adaptive filters may be required. 

Concluding Remarks 

An analysis is presented of the discontinuous 
Galerkin method for a scalar nonlinear model prob- 
lem that contains a shock. The transition from a 
stable to an unstable method as the degree of the 
method is increased from 2 to 3 is correctly pre- 
dicted by the analysis. Also, the instability asso- 
ciated with truncating the flux polynomial is also 
predicted. Insight gained from the analysis suggests 
a possible source for the instability and leads to a 
modified approximate Riemann flux that eliminates 
the instability. The elimination of the instability is 
predicted in the analysis of the modified method and 
is verified by numerical experiment using Burgers’ 
equation. Although the solutions are oscillatory in 
the neighborhood of the shock, the oscillations are 
bounded and are confined to the neighborhood of the 
shock. Furthermore, the element-averaged solution 
is monotone, which suggests that a smooth solution 
can be obtained by applying a filter as a post pro- 
cessing step. Numerical experiments with the Eu- 
ler equations showed that the modified approximate 
Riemann flux resulted in a stable scheme. Cases in 
which the shock was in motion such that eigenval- 


ues did not change sign across the shock produced 
results similar to those observed with Burgers’ equa- 
tion. In other cases, however, the oscillations were 
not confined to the shock region, and some type of 
filtering is necessary. 
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Figure 3. Test solution when shock is near right 




Figure 1. Real part, of eigenvalues of baseline dis- 
continuous Galerkin method. 




Figure 2. Real part of eigenvalues for n = 2 with Figure 4. Real part of eigenvalues of discontinuous 
flux polynomial truncated to n + 1 = 3. Galerkin method with Riemann flux given by equa- 

tion (9). 
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a) n = 2. b) n = 4. 


Figure 5. Real part of eigenvalues with averaging 
used in Riemann flux. 



o 1 u o 1 

x x 


a) Evolution of solution b) Shock at element 
with shock at element boundary: n = 0-5, t = 
boundary: n = 3. 0.4. 



c) Shock at mid-element: d) Local and element- 
77 = 0-5, t = 0.4. averaged solution with 

shock at mid-element: 
n = 3, t = 0.4. 

Figure 6. Numerical solution with modified Rie- 
mann flux for stationary shock. 



a) Local solution, b) Local and element- 

77 = 0-5, t = 0.4. averaged solution for 

?7 = 3, t = 0.4. 

Figure 7. Numerical solution with modified Rie- 
mann flux for moving shock. 



x 


b) Wco = 0.5^7. 

Figure 8. Solution to Euler equations at t = 3 with 
«co chosen such that flux splitting is essentially fixed. 
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Figure 9. Solution to Euler equations at t = 3 wit.l 
«co = ~\/l- 



X 


Figure 10. Solution to Euler equations at t = 3 wit.l 
Woo = — ^7 and n = 3 with spectral filter. 



